Introduction

In this Rmarkdown document we are going to load the MAGIC-denoised data to better visualize genes and ease with the annotation when using specific marker genes. MAGIC was developed by Smita Krishnaswamy’s lab to try to fill in the drop out reads in the spots. MAGIC is a Markov Affinity-based Graph Imputation of Cells used for denoising high-dimensional data most commonly applied to single-cell RNA sequencing data. MAGIC learns the manifold data, using the resultant graph to smooth the features and restore the structure of the data based on their k-nearest neighbors.

Libraries

library(Seurat)
library(dplyr)
library(ggplot2)
library(SPATA2)
library(UCell)
library(stringr)
library(readr)

Parameters

set.seed(123)
source(here::here("misc/paths.R"))
source(here::here("utils/bin.R"))

"{gct}/{plt_dir}" %>%
  glue::glue() %>%
  here::here() %>%
  dir.create(path = .,
             showWarnings = FALSE,
             recursive = TRUE)

"{gct}/{robj_dir}" %>%
  glue::glue() %>%
  here::here() %>%
  dir.create(path = .,
             showWarnings = FALSE,
             recursive = TRUE)

Load data

The data used in this Rmarkdown document comes from 03-clustering_integration.Rmd where the data was integrated.

merged_se <- "{anot}/{robj_dir}/integrated_spatial_annot.rds" %>%
  glue::glue() %>%
  here::here() %>%
  readRDS(file = .)

Load MAGIC data from the script MAGIC_denoising-GC.Rmd

magic_df <- "{gct}/{robj_dir}/MAGIC-mtrx.rds" %>%
  glue::glue() %>%
  here::here() %>%
  readRDS(file = .)

# create a new assay to store MAGIC information
magic_assay <- CreateAssayObject(counts = as.matrix(magic_df))

# Subset merged_se to those barcodes used
merged_se <- merged_se[, colnames(magic_df)]

# add this assay to the previously created Seurat object
merged_se[["MAGIC_Spatial"]] <- magic_assay

Analysis

Load gene list

"{gct}/GC_dict.R" %>%
  glue::glue() %>%
  here::here() %>%
  source(file = .)
gene_vec <- rownames(magic_df)
gc_vec <- unique(c(gc_dict[["HARD"]], gc_dict[["MID"]]))

Marker gene location

Look at the location where the genes of interest are expressed

Seurat::DefaultAssay(merged_se) <- "MAGIC_Spatial"

# Iterate over each image
lapply(id_sp_df$gem_id, function(i) {
  print(i)

  gene_plt <- Seurat::SpatialFeaturePlot(
    object = merged_se,
    features = gc_vec,
    alpha = c(0, 1),
    ncol = 6,
    images = i)

  # Save plot
  "{gct}/{plt_dir}/magic_GC_markers_{i}.pdf" %>%
    glue::glue() %>%
    here::here() %>%
    cowplot::save_plot(
      filename = .,
      plot = gene_plt,
      base_height = 35,
      base_width = 35)
})
## [1] "tarwe1_xott6q"
## [1] "c28w2r_7jne4i"
## [1] "esvq52_nluss5"
## [1] "p7hv1g_tjgmyj"
## [1] "gcyl7c_cec61b"
## [1] "zrt7gl_lhyyar"
## [1] "qvwc8t_2vsr67"
## [1] "exvyh1_66caqq"
## [[1]]
## [1] "/scratch/devel/melosua/phd/projects/BCLLatlas/tonsil_atlas/spatial_transcriptomics/GC-Analysis/2020-09-22/plots_2020-09-22/magic_GC_markers_tarwe1_xott6q.pdf"
## 
## [[2]]
## [1] "/scratch/devel/melosua/phd/projects/BCLLatlas/tonsil_atlas/spatial_transcriptomics/GC-Analysis/2020-09-22/plots_2020-09-22/magic_GC_markers_c28w2r_7jne4i.pdf"
## 
## [[3]]
## [1] "/scratch/devel/melosua/phd/projects/BCLLatlas/tonsil_atlas/spatial_transcriptomics/GC-Analysis/2020-09-22/plots_2020-09-22/magic_GC_markers_esvq52_nluss5.pdf"
## 
## [[4]]
## [1] "/scratch/devel/melosua/phd/projects/BCLLatlas/tonsil_atlas/spatial_transcriptomics/GC-Analysis/2020-09-22/plots_2020-09-22/magic_GC_markers_p7hv1g_tjgmyj.pdf"
## 
## [[5]]
## [1] "/scratch/devel/melosua/phd/projects/BCLLatlas/tonsil_atlas/spatial_transcriptomics/GC-Analysis/2020-09-22/plots_2020-09-22/magic_GC_markers_gcyl7c_cec61b.pdf"
## 
## [[6]]
## [1] "/scratch/devel/melosua/phd/projects/BCLLatlas/tonsil_atlas/spatial_transcriptomics/GC-Analysis/2020-09-22/plots_2020-09-22/magic_GC_markers_zrt7gl_lhyyar.pdf"
## 
## [[7]]
## [1] "/scratch/devel/melosua/phd/projects/BCLLatlas/tonsil_atlas/spatial_transcriptomics/GC-Analysis/2020-09-22/plots_2020-09-22/magic_GC_markers_qvwc8t_2vsr67.pdf"
## 
## [[8]]
## [1] "/scratch/devel/melosua/phd/projects/BCLLatlas/tonsil_atlas/spatial_transcriptomics/GC-Analysis/2020-09-22/plots_2020-09-22/magic_GC_markers_exvyh1_66caqq.pdf"

Now with the log-norm expression

Seurat::DefaultAssay(merged_se) <- "Spatial"

lapply(id_sp_df$gem_id, function(i) {
  # Iterate over each image
  gene_plt <- Seurat::SpatialFeaturePlot(
    object = merged_se,
    features = gc_vec,
    alpha = c(0, 1),
    ncol = 6,
    images = i)

  "{gct}/{plt_dir}/lognorm_GC_markers_{i}.pdf" %>%
    glue::glue() %>%
    here::here() %>%
    cowplot::save_plot(
      filename = .,
      plot = gene_plt,
      base_height = 35,
      base_width = 35)
})
## [[1]]
## [1] "/scratch/devel/melosua/phd/projects/BCLLatlas/tonsil_atlas/spatial_transcriptomics/GC-Analysis/2020-09-22/plots_2020-09-22/lognorm_GC_markers_tarwe1_xott6q.pdf"
## 
## [[2]]
## [1] "/scratch/devel/melosua/phd/projects/BCLLatlas/tonsil_atlas/spatial_transcriptomics/GC-Analysis/2020-09-22/plots_2020-09-22/lognorm_GC_markers_c28w2r_7jne4i.pdf"
## 
## [[3]]
## [1] "/scratch/devel/melosua/phd/projects/BCLLatlas/tonsil_atlas/spatial_transcriptomics/GC-Analysis/2020-09-22/plots_2020-09-22/lognorm_GC_markers_esvq52_nluss5.pdf"
## 
## [[4]]
## [1] "/scratch/devel/melosua/phd/projects/BCLLatlas/tonsil_atlas/spatial_transcriptomics/GC-Analysis/2020-09-22/plots_2020-09-22/lognorm_GC_markers_p7hv1g_tjgmyj.pdf"
## 
## [[5]]
## [1] "/scratch/devel/melosua/phd/projects/BCLLatlas/tonsil_atlas/spatial_transcriptomics/GC-Analysis/2020-09-22/plots_2020-09-22/lognorm_GC_markers_gcyl7c_cec61b.pdf"
## 
## [[6]]
## [1] "/scratch/devel/melosua/phd/projects/BCLLatlas/tonsil_atlas/spatial_transcriptomics/GC-Analysis/2020-09-22/plots_2020-09-22/lognorm_GC_markers_zrt7gl_lhyyar.pdf"
## 
## [[7]]
## [1] "/scratch/devel/melosua/phd/projects/BCLLatlas/tonsil_atlas/spatial_transcriptomics/GC-Analysis/2020-09-22/plots_2020-09-22/lognorm_GC_markers_qvwc8t_2vsr67.pdf"
## 
## [[8]]
## [1] "/scratch/devel/melosua/phd/projects/BCLLatlas/tonsil_atlas/spatial_transcriptomics/GC-Analysis/2020-09-22/plots_2020-09-22/lognorm_GC_markers_exvyh1_66caqq.pdf"
se_sub <- subset(merged_se, subset = gem_id == "esvq52_nluss5")
se_sub
## An object of class Seurat 
## 29266 features across 3079 samples within 2 assays 
## Active assay: Spatial (26846 features, 2000 variable features)
##  1 other assay present: MAGIC_Spatial
##  3 dimensional reductions calculated: pca, umap, harmony
se_sub@images <- se_sub@images[Seurat::Images(se_sub) == "esvq52_nluss5"]

Correlation matrix

Since we are working with sample esvq52_nluss5 in this example we will limit the correlation plot to this slide.

(cor_mtrx <- SCrafty::correlation_heatmap( 
  se = se_sub,
  genes = gc_vec,
  assay = "MAGIC_Spatial",
  slot = "data"))

"{gct}/{plt_dir}/magic_cor-mtrx_markers.pdf" %>%
  glue::glue() %>%
  here::here() %>%
  cowplot::save_plot(
    filename = .,
    plot = cor_mtrx,
    base_height = 15,
    base_width = 15)

  # Correlation with lognorm expression
cor_log <- SCrafty::correlation_heatmap( 
  se = se_sub,
  genes = gc_vec,
  assay = "Spatial",
  slot = "data")

"{gct}/{plt_dir}/lognorm_cor-mtrx_markers.pdf" %>%
  glue::glue() %>%
  here::here() %>%
  cowplot::save_plot(
    filename = .,
    plot = cor_log,
    base_height = 9,
    base_width = 10)

Look at them side by side

cor_mtrx + cor_log

Gene signatures

se_sub <- AddModuleScore_UCell(
    obj = se_sub,
    features = gc_dict[c("HARD", "MID", "all-targets-intersect")],
    name = '_UCell'
  )

Look at signatures

sign <- colnames(se_sub@meta.data)[stringr::str_detect(
    string = colnames(se_sub@meta.data),
    pattern = '_UCell')]

Seurat::SpatialPlot(
  object = se_sub,
  features = sign,
  alpha = c(0, 1),
  images = "esvq52_nluss5",
  ncol = 3)

Follicle focus

follicle_bc <- "{gct}/{robj_dir}/follicle_coordinates.csv" %>%
    glue::glue() %>%
    here::here() %>%
    read_csv()

interzone_bc <- "{gct}/{robj_dir}/LZ_DZ_interzone.csv" %>%
    glue::glue() %>%
    here::here() %>%
    read_csv()

se_follicle <- se_sub[, follicle_bc$barcode]
se_follicle$interzone <- colnames(se_follicle) %in% interzone_bc$barcode

Look at the signature witin the follicle

Seurat::SpatialPlot(
  object = se_follicle,
  features = sign,
  alpha = c(0, 1),
  images = "esvq52_nluss5",
  ncol = 2,
  pt.size.factor = 5) +
  Seurat::SpatialPlot(
    object = se_follicle,
    group.by = "interzone",
    images = "esvq52_nluss5",
    pt.size.factor = 5)

DE between Interzone & Follicle

Seurat::Idents(object = se_follicle) <- se_follicle@meta.data[, "interzone"]
markers_follicle <- Seurat::FindAllMarkers(
  object = se_follicle,
  assay = "Spatial",
  slot = "data",
  verbose = TRUE, 
  only.pos = TRUE)

# Look at cluster = TRUE
DT::datatable(markers_follicle, filter = "top")

Intersection with lists of interest

iz_mgs <- markers_follicle %>% filter(cluster == TRUE) %>% pull(gene)
intersect(iz_mgs, gc_dict[["HARD"]])
## character(0)
intersect(iz_mgs, gc_dict[["MID"]])
## character(0)
intersect(iz_mgs, gc_dict[["all-targets-intersect"]])
## character(0)

Session Info

sessionInfo()
## R version 4.0.1 (2020-06-06)
## Platform: x86_64-conda_cos6-linux-gnu (64-bit)
## Running under: Red Hat Enterprise Linux Server release 6.7 (Santiago)
## 
## Matrix products: default
## BLAS/LAPACK: /scratch/groups/singlecell/software/anaconda3/envs/spatial_r/lib/libopenblasp-r0.3.12.so
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C               LC_TIME=es_ES.UTF-8        LC_COLLATE=en_US.UTF-8     LC_MONETARY=es_ES.UTF-8    LC_MESSAGES=en_US.UTF-8    LC_PAPER=es_ES.UTF-8       LC_NAME=C                  LC_ADDRESS=C               LC_TELEPHONE=C             LC_MEASUREMENT=es_ES.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] readr_2.1.2        stringr_1.4.0      UCell_1.2.0        SPATA2_0.1.0       ggplot2_3.3.5      dplyr_1.0.8        SeuratObject_4.0.4 Seurat_4.1.0      
## 
## loaded via a namespace (and not attached):
##   [1] corrplot_0.92               systemfonts_1.0.4           plyr_1.8.6                  igraph_1.2.11               lazyeval_0.2.2              splines_4.0.1               crosstalk_1.2.0             listenv_0.8.0               scattermore_0.8             GenomeInfoDb_1.26.7         digest_0.6.29               htmltools_0.5.2             fansi_1.0.2                 magrittr_2.0.2              tensor_1.5                  cluster_2.1.2               ROCR_1.0-11                 limma_3.46.0                tzdb_0.2.0                  globals_0.14.0              matrixStats_0.61.0          vroom_1.5.7                 spatstat.sparse_2.1-0       colorspace_2.0-3            ggrepel_0.9.1               textshaping_0.3.6           xfun_0.30                   crayon_1.5.0                RCurl_1.98-1.6              jsonlite_1.8.0              spatstat.data_2.1-2         survival_3.3-1              zoo_1.8-9                   glue_1.6.2                  polyclip_1.10-0             gtable_0.3.0                zlibbioc_1.36.0             XVector_0.30.0              leiden_0.3.9                DelayedArray_0.16.3         future.apply_1.8.1          SingleCellExperiment_1.12.0
##  [43] BiocGenerics_0.36.1         abind_1.4-5                 scales_1.1.1                DBI_1.1.2                   spatstat.random_2.1-0       miniUI_0.1.1.1              Rcpp_1.0.8                  viridisLite_0.4.0           xtable_1.8-4                reticulate_1.24             spatstat.core_2.4-0         bit_4.0.4                   DT_0.21                     stats4_4.0.1                htmlwidgets_1.5.4           httr_1.4.2                  RColorBrewer_1.1-2          ellipsis_0.3.2              ica_1.0-2                   farver_2.1.0                pkgconfig_2.0.3             sass_0.4.0                  uwot_0.1.11                 deldir_1.0-6                here_1.0.1                  utf8_1.2.2                  ggcorrplot_0.1.3            labeling_0.4.2              tidyselect_1.1.2            rlang_1.0.2                 reshape2_1.4.4              later_1.3.0                 munsell_0.5.0               tools_4.0.1                 cli_3.2.0                   generics_0.1.2              ggridges_0.5.3              evaluate_0.15               fastmap_1.1.0               ragg_1.2.2                  yaml_2.3.5                  goftest_1.2-3              
##  [85] bit64_4.0.5                 knitr_1.37                  fitdistrplus_1.1-6          purrr_0.3.4                 RANN_2.6.1                  pbapply_1.5-0               future_1.24.0               nlme_3.1-155                mime_0.12                   SCrafty_0.1.0               compiler_4.0.1              plotly_4.10.0               png_0.1-7                   spatstat.utils_2.3-0        tibble_3.1.6                bslib_0.3.1                 stringi_1.7.6               highr_0.9                   lattice_0.20-45             Matrix_1.4-0                vctrs_0.3.8                 pillar_1.7.0                lifecycle_1.0.1             spatstat.geom_2.3-2         lmtest_0.9-39               jquerylib_0.1.4             RcppAnnoy_0.0.19            data.table_1.14.2           cowplot_1.1.1               bitops_1.0-7                irlba_2.3.5                 httpuv_1.6.5                patchwork_1.1.1             GenomicRanges_1.42.0        R6_2.5.1                    promises_1.2.0.1            KernSmooth_2.23-20          gridExtra_2.3               IRanges_2.24.1              parallelly_1.30.0           codetools_0.2-18            MASS_7.3-55                
## [127] assertthat_0.2.1            SummarizedExperiment_1.20.0 rprojroot_2.0.2             withr_2.5.0                 sctransform_0.3.3           S4Vectors_0.28.1            GenomeInfoDbData_1.2.4      mgcv_1.8-39                 parallel_4.0.1              hms_1.1.1                   grid_4.0.1                  rpart_4.1.16                tidyr_1.2.0                 rmarkdown_2.12              MatrixGenerics_1.2.1        Rtsne_0.15                  Biobase_2.50.0              shiny_1.7.1